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We present a new simulation method to calculate the free energy and the chemical po- 
tential of hard particle systems. The method rehes on the introduction of a parameter 
dependent potential to smoothly transform between the hard particle system and the 
corresponding ideal gas. We applied the method to study the phase transition behav- 
ior of monodispersed infinitely thin square platelets. First, we equilibrated the square 
platelet system for different reduced pressures with a usual isobaric Monte Carlo (MC) 
simulation and obtained a reduced pressure-chemical potential plot. Then we introduce 
the parametrized potential to interpolate the system between the ideal gas and the hard 
particles. After selecting the potential, we performed isochoric MC runs, ranging from 
the ideal gas to the hard particle limit. Through an iterative procedure, we compute 
the free energy and the chemical potential of the square platelet system by evaluat- 
ing the volume of the phase space attributed to the hard particles, and then we find 
the coexistence pressure of the system. Our method provides an intuitive approach to 
investigate the phase transitions of hard particle systems. 

KEYWORDS: multi-canonical method, isotropic-nematic phase transition, isobaric and isochoric Monte 
Carlo simulations, phase space multi-histogram method, chemical potential 

1. Introduction 

Since Onsager's pioneering work,^) it has been studied the phase transition of rod- 
shaped and disc-shaped hard particle systems using simulations with excluded volume 
effects. For example, the isotropic-nematic-smectic A phase transition was investigated 
by Frenkel et al.^^ using an isobaric MC of hard spherocy finders. Frenkel and Mulder^) 
and Veerman and Frenkel^^ observed various phase transitions of disk-shaped hard 
particles modeled as oblate elfipsoids and cut spheres respectively. The isotropic-nematic 
phase transition of monodispersed infinitely thin disks was first studied by Eppenga 
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and Prenkel,^^ and then the nematic-columbar phase transition of those was observed 
by Bates and Prenkel.^^ 

In this paper, we introduce a new simulation method which we shall call the "phase 
space multi-histogram" (PSMH) method to compute chemical potential of hard parti- 
cles. The method is a variant of the multi-canonical ensemble method/' ^-^ As an illus- 
tration of the PSMH method, we compute the variations in the chemical potentials for 
nematic and isotropic phases and detect the isotropic-nematic phase transition point 
for the square platelet system examined by Bates in 1999. 

Bates' computer simulation work^^ was motivated from an experimental result by 
van der Kooij and Lekkerkerker who synthesized suspensions of sterically stabilized 
gibbsite {Al{OH)s) platelets and observed the isotropic-nematic phase transition in a 
critical range of concentration of the experimental system. Examining the nature 
of their experimental system. Bates studied the isotropic-nematic phase transition of 
monodisperse infinitely thin platelet systems using grand canonical {fJ^VT ensemble) 
Monte Carlo (MC) simulations. It was found that the shape of a platelet is essential to 
determine its isotropic-nematic phase transition behavior. In particular, the isotropic 
and nematic coexistence reduced densities of the isotropic and nematic phases for the 
system of square platelets were in the range of 3.1-3.4 due to hysterisis. 

The organization of this paper is the foUowings. In Section 2, the PSMH method is 
explained after briefly reviewing the canonical and multi-canonical ensemble MC simu- 
lations. The method can be used to calculate the free energy of hard repulsive particles 
by interpolating from the ideal gas state. In Section 3, the PSMH method is used for the 
detection of the isotropic (1) - nematic (N) phase transition of monodispersed infinitely 
thin square platelets. Finally, the conclusion of the present study is given in Section 4. 

2. The PSMH Method 

2.1 Passing from Canonical to Multi- canonical 

Consider a non-isolated macroscopic system in contact with a heat reservior. Under 
the constraint that the number of particles N , volume V and temperature T are con- 
stant, the probability Pnvt that the system is in microstates with energy E is propor- 
tional to w{E) exp{—f3E). Here, w{E) is the number of microstates in the same energy 
E, and /3 is the inverse temperature defined by ^ = 1/(A;bT'), where ks is Bolzmann's 
constant. The ensemble defined by Pnvt is called canonical. It is known that the canon- 
ical method has the limitation of not being able to escape from a metastable state, and 
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Fig. 1. (top) Isotropic phase and (bottom) Nematic phase of the system of 120 square platelets. The 
white balls indicate the centers of mass of the square platelets. 

the multi-canonical method was introduced by Breg and Neuhaus^'®^ to remedy the sit- 
uation. The multi-canonical probability Pm has a uniform probability distribution such 
that Pm is proportional to g{E) exp{—PEfn), where g{E) is the density of states, and 
Em is the multi-canonical energy. The multi-canonical ensemble MC simulation is 
carried out on the basis of the Metropolis algorithm^^^ with the re- weighting histogram 
technique to study phase transitions.^'^) 

2.2 Introducing a Parameterized Potential 

A key ingredient of the PSMH method is a parameter dependent potential Uh to 
smoothly transform between the hard particle system and the corresponding ideal gas 
state. Recall that the classical partition function of the potential energy of a system, 
say U, is given by 



and once one knows all information about the density of states g{U), the partition 
function Z{P) is computable for any /3 multi-canonically by the re- weighting histogram 
method. Consider a hard particle system with a pairwise repulsive potential energy 
such that (pij = oo if the i- and j-th particles overlap, and (pij = if they do not. Then 




(1) 
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the total potential energy of the system is given by Uh = Ylii<j(t>ij- By Eqn. (1), the 
partition function of the potential energy Uh is written as 



= j dUhQiUh) exp{-pUh) = J dUhw(Uh)5{Uh - Uo) exp(-/3t/;,) = w{Uo), (2) 

where S is the Dirac delta function, and Uo — 0. Since the partition function Z{/3) does 
not depend on (5 any more, we may write Z{P) = Z. 

In order to compute w{0), we shall introduce a parameter dependent potential, C/^ 
defined by 

K ^ (3) 

where h is an interpolation constant, and ric is the number of pair intersections such 
that 



I 1 if the i-th and j-th particles overlap, 
otherwise. 



i<j 



In particular, for a distribution exp(— /3C/^), ric is fluctuating when h is constant, and 
we have the following relationship 

when /i — )■ oo, ric ^ (hard body model), 

(5) 

when h ^ 0, ric ^ (ideal gas). 

In reality, the parameter dependent potential C/^ behaves as the potential Uh as h goes 
to infinity 

hm U'^ = Uh, (6) 
and the partition function Z'{f3) of the parameter dependent potential U'f^ is 

Z'{l3) = j dncw{nc)ex]){-hnc). (7) 

Writing (3Uh as the product of an interpolation constant h and the number of pair 
intersections ric enables one to describe the model of an ideal gas which allows non-zero 
pair intersections. 

The PSMH method uses the re-weighting multi-histogram technique. The PSMH 
method is also similar to the multi-canonical method^' since it changes h which is 
equivalent to changing f3, but it is different in the following respect. The PSMH method 
adopts the multi-canonical method to obtain w{0); however, we do not actually perform 
the multi-canonical method for the potential Uh^oo whose partition function Z is inde- 
pendent on /3 as seen in Eqn. (2). The PSMH method is a variant of the multi-canonical 
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method^' to study phase transitions of hard particle systems with the re- weighting 
technique. 

2.3 The PSMH method: Procedures and Calculation 

First, we perform a usual isobaric MC simulation for hard particles. It corresponds 
to the limit of /i — >■ oo in Eqns. (3), (4). After equilibrating, we switch the simulation 
scheme to an isochoric one and remove the constraint of non-overlapping of hard parti- 
cles, in other words, wc perform an isochoric MC simulation with the Boltzmann weight 
factor exp(— We call the latter scheme as NVh ensemble MC simulations. In the 
NVh MC runs, the number of particle pair intersections fluctuates. Let w{nc) be the 
number of microstates such that there are Uc pair contacts in each microstate. From 
the NVh MC simulation data, we construct the apparent distribution iph{nc) of ric for a 
given h. iphinc) obeys the normalization condition '0;i(?T'c) — 1 ("Histogram"). Let 
Ch be a normalization constant and h be incremented by Ah. Then 

i/^h{nc) = Chw{nc) exp(-/3;7^), (8) 

'tph+Ah{nc) = ^^^^ V^fc(nc)exp(-A/inc). (9) 

Since ric is an integer, we call it using the letters j, k, I simply. A straightforward 
calculation gives us 

law{k)^lii^h=o{k)-lnCo,k^j,j-l,j-2,--- ,1,1 < j (10) 

when h = 0, and 

In «;(/ - 1) = In ip^il - 1) - In ^'^(0 + In w{l) + h (11) 

when h ^ 0. Here, Cq^ = w{nc). 

Assume that the global minimum value of iph{nc) is attained, and then computation 
of the free energy F and chemical potential /i of the hard particles are carried out 
via 

/3F ^ - \nw{0) +\nCo (12) 

and 



P^- N ~ N ^ N^D^ ' ^^^^ 
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where p, V, L, D, p* are a pressure, the volume of a periodic simulation box, the side 
length of the simulation box, the unit length in simulation, and a reduced pressure 

p* = D^p/3. (14) 

The estimation of Cq^ is given in Appendix , and we set the thermal de Broghe 
wavelength as A = 1 in simulation. 

3. Illustration of the PSMH method 

3.1 Model system 

Our model system is a monodispersed square platelet system. The system is com- 
posed oi N — 120 square platelets of equal area A. The unit length D in Eqn. (14) is set 
hy D = We performed isobaric (NP ensemble) MC simulations. The initial config- 
uration of the system was a dilute aligned crystal with a cubic periodic condition, and 
it was relaxed at p* = 0.1 to become an isotropic phase. We isotropically compressed it 
(i.e. the shape of the periodic simulation box was always cubic) at the reduced pressure 
ranged in p* — 0.1-2.0, where (3 is the inverse temperature. After that, we expand it at 
the reduced pressure in the same range. For compression (or expansion), an initial struc- 
ture was an equilibrium structure obtained by a previous simulation. At each pressure, 
the system was equilibrated to compute statistical quantities. The isotropic-nematic 
phase transition of the system is driven by an excluded volume repulsion (i.e. colloidal 
interactions which prohibit any overlap of platelets.) 

Following Bates' notation, we analogously define a number density p = with 
a scaled diameter a = \J^- Let p/v/ and p/yv be the number coexistence densities 
along two branches, that is, one from the nematic phase to the isotropic phase {NI 
branch) and the other from the isotropic phase to the nematic phase {IN branch) of 
the system of square platelets due to hysteresis. In simulation, the IN branch is obtained 
by compression, and the NI branch is done by expansion. 

A nematic order parameter matrix 5 is a matirx quantity which measures the 
isotropic-nematic phase transition and is given by 

1^3 1 
fc=i 

when N is the fixed number of the square platelets in the system. Cj, Cj are unit vectors 
along the lab coordinates and is a unit vector along the kth molecule's symmetric 
axis here. Since each entry {Sij) is a real number for any i,j, and (r^ • ej)(rfe • ej),Sij 
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are symmetric with respect of i and j, <S is a real symmetric matrix. Therefore, S is 
diagonahzable. The largest eigenvalue of the diagonalized S is called an order parameter 
denoted by S and its corresponding eigenvector is called a director. A model system with 
a nematic behavior shows a nematic phase with an orientational order (i.e. a director 
without any positional order.) 

3.2 Results and Discussion 

In FIG. 1, we displayed typical isotropic (S fa 0.1) and nematic (S fa 0.7) phase 
configurations of a square platelet system. Due to the finite size of the system, hysteresis 
in density was observed around p* — 1.1 as pressure was changed. (See FIG. 2.) In FIG. 
3, a plot of histograms 'tpui'i^c) for each h ranged in [0, 6] was shown. The NI branch 
at p* = 0.8 was selected initially. FIG. 4 shows the behavior of histograms when each 
h is on the interval [0,0.1]. The histograms for h = 0.1 and arc shown in dotted 
and solid curves. Wc found the histogram for /i — )■ converging to that for h = 0. 
(See FIG. 4.) The PSMH method worked in the square platelet system since multi- 
histograms could generate overlapped histograms within full width at half maximum. 
(See FIG. 5.) The ric vs. — In w{nc) plot for the NI branch at p* — 0.9 were shown in 
FIG. 6. When ric was big, — In w{nc) took a small value. — In w{nc) attained the global 
minimum at ric — 234, and we found —\nw{nc) fa 325 when ric — 0. In FIG. 7, p* vs. 
chemical potential curves for IN and NI branches was shown. The phase transition was 
observed at p* = 1.13 (^/x = 4.61). The order parameter was S f« 0.7 when the number 
coexistence densities of the system of square platelets were pjv/ ~ 3.85 and p/7v ~ 3.87 
(at p* = 1.13). This was comparable to Bates' {pni,Pin) ~ (3.1,3.4)^-* when taking 
into account of the influence of systematic errors. In colloidal suspensions of Al{OII)^ 
platelets, the number coexistence densities pm, pm and the difference between those 
densities, piN — pm were known to be very small as an experimental fact.^^ Although, in 
our simulation, the number coexistence densities were relatively high, but the difference 
was only pn^ — pm — 0.02. 

In Bates' work,^) he used a different expression: 

AttV 

Pf^Bates = I3p+ (log -^). (16) 

Our simulation result in chemical potential at the isotropic-nematic phase transition 
was Pji — 4.61, and Bates' was PjiBates — 7.3. By using Eqn. (16), we get 

(ip = PpBates - log(47ry) ^ 4. 1 (17) 
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Fig. 2. p* vs. Number density p. IN branch means a phase transition from isotropic to nematic, 
and vice versa. A weak hysteresis was observed around p* = 1.1. A vertical error bar on each data 
point indicates the standard error of the mean. For a guide to the eyes, the Unes between the data are 
drawn. 



when V is the 5x5x5 simulation box, and 

Pl^ = P^J^Bates - \0g{W) ^ 3.5 (18) 

when V is the 8x8x8 simulation box. In either case, our simulation result was 
comparable to Bates'. We did not studied the effect of the shape of the simulation box, 
but our result combined with Betes'^^ suggested that there was the size dependence 
of the free energy and chemical potential of the monodispersed infinitely thin square 
platelet system. 

While having had such a successful example of the PSMH method, the method 
has its limitation for some hard particle systems where particle insertion is rejected. 
For example, the PSMH method did not work in the sphere system to study a phase 
transition from a solid state to a liquid state. This indicates that the PSMH method 
works well in hard particle systems where particle insertion is possible. 
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Fig. 3. A plot of histograms iph{nc) for h ranged in [0,6]. The NI branch at p* = 0.8 was taken. 
When /i = 0, the maximum of tph{nc) was 0.02842 at ric = 215. 




Fig. 4. A plot of histograms t/Jhiric} for h = 0, 0.01, 0.001, 0.0001, 0.1. 
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Fig. 5. (top) A histogram ijjhinc) for h = {) and (bottom) Two histograms ijjh{nc) when h = 0. 0.1. 
The NI branch at p* = 0.8 was used, j = 319 and I = 201 were taken when h = 0. The nieghboring 
histograms overlapped within full width at half maximum. 



4. Conclusions 

We introduced the phase space multi-histogram (PSMH) method to calculate the 
free energy and the chemical potential of hard particle systems. The most essential part 
of the PSMH method is introducing the parameter dependent potential Uh to smoothly 
transform between the hard particle system {h oo) and the corresponding ideal gas 
limit (/i — )■ 0), and a basic idea of the PSMH method is shrinking the phase space of the 
ideal gas state to that of the hard particle system using the potentials Uh parameterized 
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Fig. 7. p* vs. chemical potential curves for NI and IN branches. The solid curve [NI branch) 
represents a polynomial fit to the NI data, and the dotted curve {IN branch) is a polynomial fit to 
the IN data. The phase transition was observed at p* = 1.13 = 4.61). 
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by h. 



We investigated the isotropic-nematic phase transition of the system of monodis- 
persed infinitely thin square platelets by applying the PSMH method. After equili- 
brating the square platelet system for different reduced pressures with isobaric MC 
runs, we obtained a reduced pressure-chemical potential plot. Next, we introduced the 
parametrized potential to interpolate the system from the ideal gas state to the hard 
particle system. Following that, we selected the potential, and then we performed iso- 
choric MC runs with varying from the ideal gas state to the hard particle limit. By 
the induction procedure of the PSMH, we computed the free energy and the chemical 
potential of the square platelet system by evaluating the volume of the phase space at- 
tributed to the hard particles, and then we find the coexistence pressure of the system. 
To detect an isotropic-nematic phase transition point for the square platelet system, we 
calculated the variations in chemical potentials for nematic and isotropic phases. For 
the system, the obtained number coexistence densities and the chemical potential at 
the isotropic-nematic phase transition point closely matched up with Bates' results. 
The PSMH method works well only if multi-histograms are generated so that the neigh- 
boring histograms overlap sufficiently. 

This is a pilot study, and we will continue to develop the PSMH method. In particu- 
lar, we will investigate the behavior of the system with the choice of h, the dependance 
of the size of platelets and the influence of a simulation box geometry. The next step in 
this work is to study the isotropic-biaxial nematic phase transition of monodispersed 
infinitely thin irregular platelet systems which has as of yet not been detected. 
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Appendix: Estimation of Cq^ 

In order to compute chemical potentials, we need to estimate Cq^ for a monodis- 
persed infinitely thin square platelet system. For an ideal gas. 



yN 



(A-1) 



0,ideal JSfW^N ' 
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where A = -j=L==. Here, m is the mass of a gas particle in the ideal gas. A is known 
as the thermal de Broglie wavelength. On the other hand, for platelet systems, 

'^Q,platelet " TVIA^^ ' 

Observe that Cq^ changes as V changes. 
Since 

^C'l^^teiet = Inl^^ + ln(47r)^ - InTV! - In A^^ ^ iV(ln ^ + 1), (A-3) 



we have the following expressions: 

A^A^ 



AnrV 

F = In = -^-^A^(ln ^ + 1), (A-4) 



AttV 

G = F + pV = F + r'N = -(3-'N{\n — ), (A-5) 



G , AttV 1 47r(^) 
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